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We use a high-statistic quantum Monte Carlo and Max- 
imum Entropy regularization method to compute the dy- 
namical energy correlation function (DECF) of the one- 
dimensional (ID) S = 1/2 antiferromagnetic Heisenberg 
model at finite temperatures. We also present a finite temper- 
ature analytical ansatz for the DECF which is in very good 
agreement with the numerical data in all the considered tem- 
perature range. From these results, and from a finite tem- 
perature generalisation of the mechanism proposed by Loren- 
zana and Sawatsky [Phys. Rev. Lett. 74, 1867 (1995)], we 
compute the line shape for the optical absorption spectra of 
multimagnon excitations assisted by phonons for quasi ID 
compounds. The line shape has two contributions analogous 
to the Stokes and anti-Stokes process of Raman scattering. 
Our low temperature data is in good agreement with optical 
absorption experiments of CuOs chains in Sr2Cu03. Our fi- 
nite temperature results provide a non trivial prediction on 
the dynamics of the Heisenberg model at finite temperatures 
that is easy to verify experimentally. 
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I. INTRODUCTION 

Infrared optical absorption measurementsEI on charge- 
transfer insulators such as Sr2Cu03 and the isostruc- 
tural compound Ca2Cu03 have shown that these experi- 
ments provide important insights in our understanding 
of the spin dynamics of the underlying CuOs chains. 
The observed bands are analogous to the bands found 
in cuprates with CuO planeaJlj at mid infrared energies. 

Although direct absorption of magnetic excitations is 
not possible in the usual optical transitions induced by 
dipole interaction, they are allowed by the assistance of 
phonona3 □ whose effect is to effectively lower the symrae.- 
try of the lattice giving rise to a finite dipole momenttj. 
This effective dipole moment for processes involving one 
phonon and spin excitations has been obtained pertur- 
batively starting from the three-band Hubbard model in 
the presence of slowly varying electric and phonon fields 
(see Ref. hereafter LS) or using exact diagonalization 
of small clusteralil. 

An effective low-energy description of the spin de- 
grees of freedom of the CuOa chain is given by the one- 



dimensional (ID) S 
Hamiltonian (AFH). 



1/2 antiferromagnetic Heisenberg 



(1) 



Although the energy spectrum of this many-bodv_fluan- 
tum problem has been solved a long time agal3, the 
computation of its dynamical properties remains a chal- 
lenge. So far most of the numerical results on dynami- 
cal properties come from exact diagonalization of small 
systems, and more recently from quantum Monte Carlo 
(QMC) simulations combined with the Maximum En- 
tropy (MaxEnt) analytic continuation procedure of noisy 
data. This combination has been used to cal c . uLat e the 
dynaynjcal spin correlation function of the IDllaEJ and 
2Dt§E3 Heisenberg model and correctly describes the 
inelastic pfijitron scattering of CuCl2.2N(C5H5)E3 and 
La2Cu04Bn. 

The estimated-iatrachain exchange energy of Sr2Cu03 
is J '--^ 0.26eytH3llJ. For an electric field parallel to the 
CuOa chain the absorption band at low temperatures 
(T ~ 32K) is wide, starts at the relevant phonon fre- 
quency (wo = O.OTeV), extends up to w = + ttJ 
0.9eV (the maximum energy of two spinous dispersion) 
and has a singularity at w = -f 7rJ/2 ~ 0.48eV (the 
maximum of the des Cloizeaux-Pearson dispersion) . For 
an electric field perpendicular to the chain axis, the spec- 
trum is weak and structureless. This kind of response is 
typical of ID systems, p 

Experimental featuresfil of the optical absorption spec- 
trum a(a;) of these materials may be analysed applying 
LS idea of phonon-assisted multimagnon optical absorp- 
tion with a T = analytical ansatz for the dynamical 
energy correlation function (DECF)eI. Numerical studies 
show that at T = the dynamical structure factor and 
the DECF of the AFH model have a significant structure 
in the uj — q plane between u)i{g) = 7rJ| sin((7)|/2 (the 
des Cloizeaux-Pearson dispersion relation) and cj2('7) = 
■kJ\ sm[q)\, the maximum of the two-spinon continuum. 
At zero temperature the LS mechanism gives structure 
at the energy of the relevant phonon plus the energy of 
the magnetic excitations and the computed spectrum is 
in excellent agreement with the experimental oneu at low 
temperatures. 

Here we present numerical results for the finite tem- 
perature DECF. For systems up to Ng — 32, we perform 



1 



a high-statistic quantum Monte Carlo simulation and ob- 
tain the spectrum by the MaxEnt analytic continuation 
regularization technique of noisy data. The Monte Carlo 
technique is a non-perturbative approach and so useful 
to check the validity of the different approximations in- 
volved on the analytical results; the considered system 
sizes are large enough to safely extrapolate to the ther- 



modynamic limit. The technique is presented on Sec. Ill 



One of our goals is to study finite temperature effects 
on the optical absorption spectra. Although the temper- 
ature dependence of the speejtra has been studied exper- 
imentally in 2D compoundMi3 to the best of our knowl- 
edge no ID studies have been reported yet. In Sec. |^ we 
discuss LS theory at finite temperatures. p 

We also check the validity of the T = aasatzS for 
the DECF and, based on bosonization resultslla, we gen- 
eralise this ansatz to finite temperatures (Sec. The 
analytical ansatz is in very good agreement with the nu- 
merical data in all the considered temperature range. Fi- 
nally, our results for the optical absorption spectra (Sec. 
^) are in excellent agreement with available low temper- 
ature experimental (Ref. ^ data and provide a prediction 
on the finite temperature dynamics of a ID Heisenberg 
system that can be tested experimentally. 



II. PHONON-ASSISTED IR ABSORPTION OF 
MAGNONS AT FINITE TEMPERATURES 

A. Model 

In this section we generalise the T = theory of 
LS to finite temperatures. We consider an effective 
Hamiltonian describing the coupling of light with onc- 
phonon-multimagnon excitations: H — Hq + Hi, where 
Hi = —PE, P (E) is the magnitude of the dipole mo- 
ment operator (electric field) in the chain direction. The 
unperturbed Hamiltonian is: 



Hq — H„ 



H 



ph- 



(2) 



here the magnetic part is the Heisenberg model Eq. (|l|) 
and for the phonon Hamiltonian we take a singleJpranch: 
the Cu-0 stretching mode in the chain directiontl 



(3) 



where uj„ is the phonon frequency. The dipole moment is 



(4) 



where Aq = (a^ + aL^), Bq is the Fourier transform of 
Bi = SiSi+i. The strength of the coupling of light with 
these excitations is given by: 



(5) 



here M is the oxygen mass and qa a material dependent 
effective charge which is in the range 0.1 ^ 0.4eErS. 

In the definition of the dipole moment operator we 
have neglected terms like conventional phonon absorp- 
tion (P = X'Aq) which are outside the scope of this work. 



B. Optical absorption 



From general considerations the optical absorption 
spectrum is given by 



2uj 



(6) 



where e = ei(w) + *e2('^) is the complex-dielectric func- 
tion and c is the speed of light. Assuming weak absorp- 
tion. El » £2, Eq. (13) can be rewritten as. 



47r 



(7) 



where (t{uj) is the frequency dependent conductivity 
which can be written in terms of the Fourier transform 
of the real time dipole-dipole retarded Green function. 
Using the Lehman representation and the fluctuation- 
dissipation theorem one can write for the real part, 

^^TTu^O^^e^ ^ e-P^«\{N\P\M)\'Siu; + ej,-eM) 



ZV 



(8) 



here Hq\M) — eM\M), Z is the partition function, 
the inverse temperature and V the volume. In Eq. (g) 
weak magnon-phonon interactions were assumed, i.e in 
the absence of electric fields any term coupling magnons 
with phonons was neglected. Under this assumption the 
eigenstates of the system are products of phonon states 
times magnetic states. 



|Af) = |m)|M) 



(9) 



where \^) (|m)) labels an eigenstate of i?mag (-ffph) with 
energy E^ (wm). (From now on we use ^,,v to label 
magnetic eigenstates and m, n for phonon eigenstates.) 
The optical conductivity can be written as: 



V ZirisigZph 



J2 ^-0(^M)y2^ (10) 



\{n\Aq\7n)\^\{ll\Bq\u)\^5{LO + Un - UJn 



E,, — E,i 



where Zph = J2n e and Z^ag ^ 



Using 



the fact that Aq creates or destroys a phonon of momen- 
tum q, we can evaluate the phonon part and we get for 
the optical conductivity, 
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V " 



the method divides the < r < /? interval into L parts, 
each one of width At — (3/L and breaks the ID Hamil- 



y + nq)Jq{Lo - LOq) + 71, J,(w + ujq)] (11) ^^^ian ff^ag m two commuting parts ffi and H2, each 



Ug is the Bose occupation number for g-momentum 
phonons. All the magnetic information is stored in the 
function, 



^mag 



(12) 



The first term in the brackets in Eq. ( |11D corresponds 
to process in which the incoming photon creates a mag- 
netic excitation and a phonon; the second term corre- 
sponds to the destruction of a phonon from the bath with 
creation of a magnetic excitation. The first term does not 
vanish at zero temperature and it is the one considered by 
LS. It gives structure at the energy of the phonons plus 
the energy of the magnetic excitations. The last term is 
new, vanish at zero temperature and gives structure at 
the energy of the magnetic excitation minus the energy 
of the phonon. Notice the similarity with the Stokes and 
anti-Stokes processes in Raman scattering. 

Jq can be evaluated from the retarded Green function 
using the fluctuation-dissipation theorem: 



Jq{uj) = 



— Im Gqioj) 



7r(l 



(13) 



where Gq{u>) is the Fourier transform of the real time 
retarded Green function 

Gq{t) = -l9{t) < [Biq, t),Bi-q, 0)] > . (14) 

We also define 

Aqiu;) = (1 - e-'^")Jg(cc;) = --lmGq{Lo) (15) 

TT 

which is an odd function of lu and is evaluated in the next 



section using MaxEnt and in Sec. IV using an analytical 
ansatz. 



III. QUANTUM MONTE CARLO AND 
MAXIMUM-ENTROPY ANALYTIC 
CONTINUATION 

In order to compute the magnetic response we pQi^ 
form a World Line quantum Monte Carlo simulation.tj 
The algorithm is based on a path integral approach. To 
evaluate the average of physical observables. 



(0) = 



tr(Oe 



) 



^mag 



one corresponding to a different sublattice, even and odd 
respectively. Introducing at each one of the 2L interme- 
diate times a complete set-of Sz eigenstates and using 
the Suzuki- Trotter formulae^ 



-AtH„ 



-AtH,-AtH2 



1 + 0(At2)], 



one can evaluate the resulting matrix elements and in- 
sure, making At — *■ 0, that the error introduced by the 
Suzuki- Trotter approximation is smaller than the statisti- 
cal error. The updating procedure has two different kinds 
of moves: local changes conserve the total spin projectioii 
Sz and global moves allow magnetisation fluctuations.Ed 

We implemented a multispin code. This allows an effi- 
cient integer algebra manipulation of one-bit spins packed 
into 32-bits words and achieves a high speed. For differ- 
ent values oi /3 — 1, 5, 10 (in units of J) we fix At = 0.125; 
given this value of At the acceptance rates of local and 
global moves are always less than 10%, this is fixed by 
the algorithm and can not be changed. For /3 = 10 we 
performed ~ 100 independent runs, 10^ iterations each, 
measuring every 1000. It was necessary to evaluate the 
spin autocorrelation function in order to insure that mea- 
surements were not correlated. 

Quantum Monte Carlo simulations allow us to com- 
pute imaginary-time t correlation functions. Besides the 
DECF results presented in this work, as a check of our 
code, we have also computed static properties, the spin 
dynamic structure factor and for small systems, up to 
Ns = 16 sites, tested against exact diagonalizationca 
obtaining the optical absorption spectra from the spec- 
tral representation of the DECF. All the results shown 
here correspond to a system of size Ng = 32. For 
low temperatures, /3 = 10 and Ng = 32 we monitored 
(E)^ - {E^) = ^ Ei=o,...,iv,-i < B^B,+i >: in aU cases 
the differences were less than 5 x 10"''. 

The imaginary time energy correlation function is pe- 
riodic in T with period f3 and is given by 



Gq{T) - {Bq{T)B^q{0)) 



(16) 



where Bq{T) = e^^BqC^''^ and < t < /3. 

For the analytic continuation problem of calculating 
the spectral density Aq(a;) associated with the imaginary- 
time Green's function Gq (t) one has to invert the integral 
equation, 



Gq(r) = 



-t-00 



dujK{T, U!)Aq{uj) 



where K is the Kernel for bosons 



1 



-f3LU 



(17) 



(18) 



It is well known that Fredholm integral equations of the 
first kind are classical examples of ill-posed problems. Af- 
ter discretization of the frequency interval one ends up 
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with a discrete ill-posed problem. Because QMC simula- 
tions provide a noisy and incomplete set of data Gq(r„), 
T„ = uAt the direct inversion is ill-conditioned and the 
Ag(a;) spectrum cannot be uniquely determined. Here 
the bar over Gq{Tn) indicates a Monte Carlo average as 
explained below. 

The MaxEnt approach is a regularization method for 
this kind of ill-posed problems which provides an unique 
solution ccmpatible with Bayesian statistic and prior 
knowledge.ES The best image compatible with the data 
is obtained maximising 



AS* 
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by solving VQ — 0. Here A is a regularization parameter 
and the misfit statistic which not only includes the 
statistical errors of the data but also the correlations be- 
tween different imaginary-time data points. The misfit 
statistic is given by 

X^ = V[G,(n) ~ G,{n)]C-MG,{T,) - G,{t,)] (19) 



-Prob. 
- 90"/o 



Kolmogorov-Smimov ^ \ 
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where G is the covariance matrix that takes into account 
the imaginary-time correlation of the data. S is the 
Shannon-Jaynes entropy of the spectrum defined with 
respect to the initial guess model m{u!): 



S 



+ 00 



du;[Aq{u;) ~ m{uj) — Aq{Lu)ln{Aq{uj) /m{Lu))] 

(20) 



m{uj) is the best guess in absence of prior knowledge on 
the spectrum. We took a flat model. 

In practice, functional maximisation is reduced to 
search in a multidimensional parameter space by a Simp- 
son evaluation of the integrals. We performed ar-i^ariable 
w-discretization with 100-200 frequency points.E3 

We binned the data Gq until we obtained a Gaussian 
distribution that satisfies several criteria, Kolmogorofc 
Smirnov, skewness, Kurtosis and standard deviation.E^I 
Fig. 1^ presents an example: starting from ~ 100000 
QMC-data the four gaussianity tests are satisfied group- 
ing the data on 100 — 1000 bins, each one an aver- 
age of 1000 — 100 QMC-data respectively. We consider 
N — 1000 bins. After this binning the imaginary time 

Monte Carlo averages Gq have relative errors less than 
10-3 

A zero temperature sum rulei was used to check the 
low temperature data and was found to be satisfied 
within a 2% error. 



FIG. 1. The 4 difFerent gaussianity tests that we per- 
formed on our data. The independent variable is A'" the num- 
ber of bins. For Kolmogorov-Smirnov we plot the probability 
of the data being Gaussian. We can see that up to A'' = 1000 
bins there is a high probability of the data being Gaussian. 
For Skewness and Kurtosis crossing over the full line indicates 
departure from Gaussianity. From the standard deviation (a) 
test we can see that under N = 100 bins the data does not 
behave as Gaussian. 



IV. ANSATZ FOR THE DYNAMIC ENERGY 
CORRELATION FUNCTION AT FINITE 
TEMPERATURES 

We present here a finite temperatures generalisation of 
the T = ansatz (Ref. which agrees very well with 
the numerical results. 

The T = ansatz can be written in terms of G^^ [u), 
the retarded Green function of the XY model. 



^ T nXY't N _ ^[^ - ^i(g)]^[^2(g) - iv]yjujl{q) - 

(21) 

here the prime indicates that the energy scale has been 
changed according to Jxy Jtt/2. With this notation 
the zero temperature ansatz reads 



ImG,(c^) = ImGf^ iu)[A^JU^ + B J fq{Lo)] (22) 
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where A and B are constants of values A= 
and we defined 



1 



^Cj2 - Cj2(g) 



2.4, B= 0.6 



(23) 



The term - ujf{q) in Eqs. (||),(|2|) is chosen to 

reproduce the behaviour of bosonization at low energies. 
For details see Ref. ^. 

To generalise this ansatz to finite temperatures we 
simply replace {uj) and fq{u!) by their finite tem- 



perature versions. At finite temperature it is trivial to 
compute G^^ iw) 
Fermion modeta, 



compute G^^ iw) by mapping the XY model to a free 



ImG 



XY' 



cosh(%^)+cosh(^^g|g^ 



(24) 



On the other hand the low energy behaviour predicted 
by bosonization bccomeal3: 

f,i.) ^ -2MmP^^^]W^^^]) (25) 



2tt 



2tt 



where we have defined 



h{y) 



1 ^(i + i^y) 



8^r(| + iiy) 



(26) 



Eqs. (|22|),(P4|)-(p6D define our finite temperature ansatz. 
The spectral density A is computed from the Green func- 
tion through expression Eq. (|l5|). 




FIG. 2. The spectral density A as a function of q and uj for 
/3 = 10 (J = 1). The dots are the MaxEnt data for Ns = 32 
and the lines the finite temperature analytic ansatz. 



In Fig. g we compare the finite temperature ansatz 
with the Maximum Entropy data for different momenta 
and temperatures. We see that both approaches give 
quite similar results. In particular the width of the peak 
close to q = TT agrees remarkably well (See Fig. m. 



k<-/?=10 



A 




FIG. 3. A as function of uj for g = vr and /3 = 5, 10 (J = 1). 
Dots and lines as in Fig. 

For q ^ the ansatz gives a much narrower peak than 
the MaxEnt data, however this disagreement is of no im- 
portance for our purpose since this part of the spectrum 
gives a negligible contribution to the IR line shape. 



V. INFRARED ABSORPTION SPECTRA 

The optical spectra for phonon-assisted multimagnon 
absorption is given by Eqs. (fzl) , (|ll|) . For the materi- 
als under consideration the momentum dependence of 
the phonons can be neglected and we can take Einstein 



phonons: 



UJQ, Uq = no- 



In this case it is useful to 



define the following function. 



— J2,in\q/2)lmGq{L0), 



(27) 



In Fig. I we show this function for different values of 
the inverse temperature. From this figure the absorption 
coefficient for other materials can be computed without 
the cumbersome integration of the ansatz. The absorp- 
tion coefficient for a specific material is, 



Q!0<^(1 



-/3^Nr (l +"0)-^(^ -^O) I ngljuJ+UJo) 



where 



^Ns^'qj 
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FIG. 4. /(t^) for different values of /3 in units such that 
J = 1. 

In Fig. ^ we show our results for the midinfrared ab- 
sorption spectra obtained from the analytical ansatz and 
compare with experimental measurements on Sr2Cu03 at 
T = 32K. As a check we can see that our low tempera- 
ture results are in excellent agreement with experimental 
data. 




FIG. 5. We show a as a function of uj for Sr2Cu03 
(J = 0.246eV, ujo = O.OSeV) for (from bottom to top) 
/3/J = 5, 10, 30, 100, using the finite temperature ansatz. The 
curves are shifted by 0.1 for clarity. The thick line is the 
experimental data after Ref. m 



VI. SUMMARY AND DISCUSSION 



numerical data in all the temperature range. Since these 
two approaches are independent the agreement gives us 
confidence on the accuracy of both results. In particular 
for q = TT, the agreement of the width of the spectral func- 
tion is a quite strong test for the numerical data of finite 
systems because in the ansatz that width is completely 
determined by bosonization (which is an asintotically ex- 
act approach). 

From the ansatz we obtain the optical absorption 
which, at low temperatures, is in excellent agreement 
with available CuOs measurements. This is expected 
since the finite temperature ansatz converges to the zero 
temperature ansatz as the temperature is lowered and for 
the latter good agreement has already been foundcl. 

Naively one can expect that finite temperature effects 
will affect the line shape at energies smaller than the 
temperature. As Fig. ^ show this expectation is wrong. 
Dramatic changes occur close to the singularity at en- 
ergy luq + TT J/2. This can be understood by considering 
the mapping pOf the problem to a half-filed interacting- 
fermion bandcj. The particle hole transitions that have 
the required energy to contribute to the singularity are 
excitations of fermions from just below the Fermi level to 
the top of the band and particle hole transitions from the 
bottom of the band to just above the Fermi level. Since 
both process involve the fermions at the Fermi level they 
get strongly affected by the smoothing of the Fermi distri- 
bution function around the Fermi energy. This explains 
the rounding of the singularity with temperature. Of 
course the low energy part of the spectrum also involves 
excitations close to the Fermi level and gets affected as 
can be seen from Fig. IJ. 

Our finite temperature line shape is a non trivial pre- 
diction on the dynamics of the ID Heisenberg model at 
finite temperatures. We hope that this result will stimu- 
late further experimental work in this direction. 
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We have used quantum Monte Carlo and MaxEnt 
methods to calculate the dynamical energy correlation 
function, Aq{uj), for S= 1/2 antiferromagnetic Heisen- 
berg chains at finite temperatures. 

Based on bosonization at finite temperatures we also 
present a finite temperature ansatz for the spectral den- 
sity of the DECF which is in very good agreement with 
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